/*
 * DenseCholesky.cs
 * 
 * Copyright (c) 2006, dnAnalytics. All rights reserved.
*/
using System;
using dnAnalytics.Math;
using dnAnalytics.Resources;
using dnAnalytics.Exceptions;

namespace dnAnalytics.LinearAlgebra.Decomposition
{

    /// <summary>
    /// The Compute method locks the object and call InternalCompute. InternalCompute will 
    /// decompose the matrix, overwriting the current matrix with the decomposition, and sets the
    /// the value of isPositiveDefinite.
    /// </summary>
    internal sealed class DenseCholesky : Algorithm, ICholesky
    {
        #region Fields

        private bool isPositiveDefinite;
        private DenseMatrix matrix;

        #endregion Fields

        #region Constructor
        /// <summary>
        /// Constructs a DenseCholesky object.
        /// </summary>
        /// <param name="matrix">The matrix to decompose.</param>
        public DenseCholesky(Matrix matrix)
        {
            if (matrix == null)
            {
                throw new ArgumentNullException("matrix", Strings.NullParameterException);
            }

            if (matrix.Rows != matrix.Columns)
            {
                throw new NotSquareMatrixException(Strings.MustBeSquare);
            }

            this.matrix = new DenseMatrix(matrix);
        }

        #endregion Constructor

        #region Public Methods

        ///<summary>Return a value indicating whether the matrix is positive definite.</summary>
        ///<returns><b>true</b> if the matrix is positive definite; otherwise, <b>false</b>.</returns>
        public bool IsPositiveDefinite()
        {
            Compute();
            return isPositiveDefinite;
        }

        ///<summary>Returns the Cholesky factored matrix (lower triangular form).</summary>
        ///<returns>The lower triangular Cholesky factored matrix.</returns>
        ///<exception cref="NotPositiveDefiniteException">If the matrix is not positive definite.</exception>
        public Matrix Factor()
        {
            return matrix.Clone();
        }

        /// <summary>
        /// Copies the Cholesky factored matrix (lower triangular form) into the result <see cref="Matrix"/>.
        /// </summary>
        /// <param name="result">The <see cref="Matrix"/> store factored matrix into.</param>
        /// <exception cref="ArgumentNullException">If <paramref name="result"/> is <b>null</b>.</exception>
        /// <exception cref="NotSquareMatrixException">If <paramref name="result"/> is not a square matrix.</exception>
        /// <exception cref="NotPositiveDefiniteException">If the matrix is not positive definite.</exception>
        public void Factor(Matrix result)
        {
            matrix.CopyTo(result);
        }

        /// <summary>Calculates the determinant of the matrix.</summary>
        /// <returns>The determinant of the matrix.</returns>
        /// <exception cref="NotPositiveDefiniteException">If the matrix is not positive definite.</exception>
        public double Determinant()
        {
            double ret = 1.0;
            int rows = matrix.Rows;
            for (int j = 0; j < rows; j++)
            {
                double value = matrix.data[j * rows + j];
                ret *= (value * value);
            }
            return ret;
        }

        ///<summary>Solves a system on linear equations, AX=B, where A is the factored matrix.</summary>
        ///<param name="rightHandSide">RHS side of the system.</param>
        ///<returns>the solution matrix, X.</returns>  
        ///<exception cref="ArgumentNullException">B is null.</exception>
        ///<exception cref="NotPositiveDefiniteException">A is not positive definite.</exception>
        ///<exception cref="ArgumentException">The number of rows of A and B must be the same.</exception>
        public Matrix Solve(Matrix rightHandSide)
        {
            Matrix result = MatrixBuilder.CreateMatrix(matrix.Rows, MatrixType.Dense);
            Solve(rightHandSide, result);
            return result;
        }

        /// <summary>
        ///  Solves a system of linear equations, <b>AX = B</b>.
        /// </summary>
        /// <param name="rightHandSide">The right hand side <see cref="Matrix"/>, <b>B</b>.</param>
        /// <param name="result">The left hand side <see cref="Vector"/>, <b>x</b>.</param>
        public void Solve(Matrix rightHandSide, Matrix result)
        {
            DenseMatrix rhs = new DenseMatrix(rightHandSide);
            matrix.lapack.CholeskySolve(matrix.Rows, rightHandSide.Columns, matrix.data, rhs.data);
            rhs.CopyTo(result);
        }

        /// <summary>
        /// Solves a system of linear equations, <b>Ax = b</b>.
        /// </summary>
        /// <param name="rightHandSide">The right hand side vector, <b>b</b>.</param>
        /// <returns>The left hand side <see cref="Vector"/>, <b>x</b>.</returns>
        public Vector Solve(Vector rightHandSide)
        {
            Vector result = VectorBuilder.CreateVector(matrix.Rows, VectorType.Dense);
            Solve(rightHandSide, result);
            return result;
        }

        /// <summary>
        /// Solves a system of linear equations, <b>Ax = b</b>.
        /// </summary>
        /// <param name="rightHandSide">The right hand side <see cref="Vector"/>, <b>b</b>.</param>
        /// <param name="result">The left hand side <see cref="Vector"/>, <b>x</b>.</param>
        public void Solve(Vector rightHandSide, Vector result)
        {
            double[] rhs = rightHandSide.ToArray();
            matrix.lapack.CholeskySolve(matrix.Rows, 1, matrix.data, rhs);
            result.SetValues(rhs);
        }

        #endregion  Public Methods


        #region Protected Members
        protected override void InternalCompute()
        {
            int ret = matrix.lapack.CholeskyFactor(matrix.Rows, matrix.data);
            if (ret == 0)
            {
                isPositiveDefinite = true;
            }
        }

        #endregion Protected Members
    }
}
